1 OVERVIEW

Question this study addresses:

NOTE: THE ANALYSES CONTAINED HEREIN WILL BE FOR A MANUSCRIPT BASED ON APS 2017.
ADDITIONAL ANALYSES AND RESULTS FOR THE POSTER WILL BE IN A DIFFERENT FILE.

FIX THE TEXT BELOW

1.1 ABSTRACT

1.1.1 Title

Significant differences in sleep and depression among chronic pain groups: Implications for patient education

1.1.2 Purpose

2 START

3 LIBRARIES

# library(multcomp)
# library(knitr)
# library(readr)
library(formatR)
library(caret)
library(kableExtra)
library(psych)
library(corrplot)
library(janitor)
library(ggpubr)
library(ggsignif)
library(plyr)
library(reshape2)
library(tidyverse)
library(reprex)
library(here)
here()
dr_here(show_reason = FALSE)

4 DATA

Import data from GitHub

#       IMPORT DATA FROM GITHUB
data <- read.csv(url('https://raw.githubusercontent.com/BrainStormCenter/Manuscripts/master/Mood_DMN_2gps/behavioral_data_13042020.csv'), header = TRUE)
data <- as_tibble(data)

# 'https://raw.githubusercontent.com/BrainStormCenter/Manuscripts/Mood_DMN_2gps/behavioral_data_13042020.csv'), header = TRUE)
#https://raw.githubusercontent.com/BrainStormCenter/Manuscripts/master/Mood_DMN_2gps/behavioral_data_13042020.csv

# Mood_DMN_2gps/behavioral_data_13042020.csv

Import Subs list from GitHub

subs <- read.csv(url('https://raw.githubusercontent.com/BrainStormCenter/Manuscripts/master/Mood_DMN_2gps/SubjectList_2020_12_11_v1.csv'), header = TRUE)
subs <- as_tibble(subs)
subs <- subs %>% dplyr::select("ID") # has imaging data

4.1 Subset

#   CREATE DATA SUBSET WITH ONLY SUBJECTS WITH IMAGING DATA
dat0 <- merge(data, subs)
dat0 <- as_tibble(dat0)
# CLEANING AND PRESERVING ORIGINAL DATA

data %>% count(group)
# A tibble: 5 x 2
  group     n
  <int> <int>
1     0     1
2     1    53
3     2    80
4     3    52
5    NA     1
dat1 <- filter(dat0, group != 0, group != "NA", dat0$bdi1 != "NA")
dat1$group.factor = factor(dat1$group, levels = c("1", "2", "3"))
levels(dat1$group.factor) = c("HC", "CLBP", "FM", na.rm = TRUE)
# names(dat1)[names(dat1) == 'gp.fct'] <- 'group.factor' REPLACE HC McGILL TOTAL
# 'NA' WITH 0
dat1$mcgill_total <- coalesce(dat1$mcgill_total, 0L)

## CLEAN ESS BY REMOVING 'NA' ##
dat1 <- filter(dat1, ess_total != "NA")

## CHANGE VARIABLE TYPE ##
dat1$ess_total <- as.numeric(dat1$ess_total)
dat1$bdi_total <- as.numeric(dat1$bdi_total)
dat1$mcgill_total <- as.numeric(dat1$mcgill_total)

## ISI ## The Insomnia Severity Index score 'isi_total' is already computed


## REDUCE GROUPS TO JUST HC AND CLBP ## dat1 <- filter(dat1, group.factor != 'FM')

## FINAL SUBJECT COUNT ##
dat1 %>% count(group.factor)
# A tibble: 3 x 2
  group.factor     n
  <fct>        <int>
1 HC              35
2 CLBP            63
3 FM              25
## END ##

4.2 PSQI

  • The global PSQI score is then calculated by totaling the seven component scores, providing an overall score ranging from 0 to 21, where lower scores denote a healthier sleep quality.
  • A PSQI Global score ≥ 5 is the clinical cutoff score (i.e., sleep problems)
#       PSQI DATA
# gather psqi variables into separate data frame
dat.psqi <- dat1[,c(1,425,769,
                      grep("pittsburgh", colnames(dat1)),
                      grep("psqi", colnames(dat1)))]
#grep("timestamp", colnames(dat.raw))

# colnames(dat.psqi)
# label(dat.psqi)
# view(dat.psqi)

#   REORGANIZE VARIABLES
dat.psqi <- dat.psqi %>% select(1:7,43,everything()) 
dat.psqi <- dat.psqi %>% select(1:11,44,everything()) 
dat.psqi <- dat.psqi %>% select(1:14,45,everything()) 
dat.psqi <- dat.psqi %>% select(1:16,46,everything()) 
dat.psqi <- dat.psqi %>% select(1:18,47,everything())
dat.psqi <- dat.psqi %>% select(1:20,48,everything())
dat.psqi <- dat.psqi %>% select(1:22,49,everything())
dat.psqi <- dat.psqi %>% select(1:24,50,everything())
dat.psqi <- dat.psqi %>% select(1:26,51,everything())
dat.psqi <- dat.psqi %>% select(1:28,52,everything())
dat.psqi <- dat.psqi %>% select(1:30,53,everything())
dat.psqi <- dat.psqi %>% select(1:32,54,everything())
dat.psqi <- dat.psqi %>% select(1:35,55,everything())
dat.psqi <- dat.psqi %>% select(1:37,56,everything())
dat.psqi <- dat.psqi %>% select(1:39,57,everything())
dat.psqi <- dat.psqi %>% select(1:41,psqi_8b.factor,
                                psqi_9, psqi_9.factor,
                                everything())
dat.psqi <- dat.psqi %>% select(1:45,
                                psqi_10.factor,
                                psqi_10a, psqi_10a.factor,
                                psqi_10b, psqi_10b.factor,
                                psqi_10c, psqi_10c.factor,
                                psqi_10d, psqi_10d.factor,
                                psqi_10e, 
                                psqi_10e1, psqi_10e1.factor,
                                everything())

dat.psqi <- dat.psqi %>% select(1:2,group.factor,everything())


####            chr to num, eval=TRUE, echo=FALSE}          ####
dat.psqi$psqi_2 <- as.numeric(dat.psqi$psqi_2)
dat.psqi$psqi_4 <- as.numeric(dat.psqi$psqi_4)
dat.psqi$psqi_5a <- as.numeric(dat.psqi$psqi_5a)
dat.psqi$psqi_5b <- as.numeric(dat.psqi$psqi_5b)
dat.psqi$psqi_5c <- as.numeric(dat.psqi$psqi_5c)
dat.psqi$psqi_5d <- as.numeric(dat.psqi$psqi_5d)
dat.psqi$psqi_5e <- as.numeric(dat.psqi$psqi_5e)
dat.psqi$psqi_5f <- as.numeric(dat.psqi$psqi_5f)
dat.psqi$psqi_5g <- as.numeric(dat.psqi$psqi_5g)
dat.psqi$psqi_5h <- as.numeric(dat.psqi$psqi_5h)
dat.psqi$psqi_5i <- as.numeric(dat.psqi$psqi_5i)
dat.psqi$psqi_5othera <- as.numeric(dat.psqi$psqi_5othera)
dat.psqi$psqi_7 <- as.numeric(dat.psqi$psqi_7)
dat.psqi$psqi_8a <- as.numeric(dat.psqi$psqi_8a)
dat.psqi$psqi_9 <- as.numeric(dat.psqi$psqi_9)
dat.psqi$psqi_component1 <- as.numeric(dat.psqi$psqi_component1)


#       REPLACE "NA" WITH 0
dat.psqi$psqi_5othera <- dat.psqi$psqi_5othera %>% replace_na(0)


####            PSQI TIB SANITY CHECK           ####
# max(dat1$TIB)

# colnames(dat1)
# da# z <- data[which(dat.psqi$TIB > 19, ) ]
#       Estimate time in bed (TIB)

dat.psqi$bedtime  <- dat.psqi$psqi_1
dat.psqi$bt.sleep <- dat.psqi$psqi_1a.factor
dat.psqi$waketime <- dat.psqi$psqi_3
dat.psqi$wt.wake <- dat.psqi$psqi_3a.factor

#       FOR FORUM HELP VERSION #1
#       GOING TO SLEEP AND WAKING UP ON THE SAME DAY RESULTS IN 24+ HOURS
dat.psqi <- dat.psqi %>%
    mutate_all(.funs = as.character) %>%
    mutate_at(.vars = c("bt.sleep", "wt.wake"),
              .funs = ~ gsub(pattern = ".",
                           replacement = "",
                           x = .x,
                           fixed = TRUE)) %>%
    mutate(Sleep = as.POSIXct(x = paste(bedtime, bt.sleep),
                              format = "%I:%M %p"),
           Wake = as.POSIXct(x = paste(waketime, wt.wake),
                          format = "%I:%M %p") +
            as.difftime(tim = 1, units = "days"),
           `Time in bed` = Wake - Sleep) # TIB version 1

#       CRUDE FIX FOR THE 24+ HOUR ISSUE
dat.psqi$TIB  <- as.numeric(dat.psqi$`Time in bed`)   # TIB version 2

# round(dat.psqi$TIB  <- if_else(dat.psqi$TIB < 24,
#                              dat.psqi$TIB,
#                              dat.psqi$TIB - 24), 2)

# 
# dat.psqi$TIB2  <- if_else(dat.psqi$TIB < 24,
#                              dat.psqi$TIB,
#                              dat.psqi$TIB - 24) # TIB version 3
# 
# 


# ###           Second version from forum       ###
# dat.psqi <- dat.psqi %>%
#   mutate(Sleep = if_else(
#       bt.sleep == "PM", true = as.POSIXct(x = paste(bedtime, bt.sleep), 
#                                           format = "%I:%M %p"), # time in current date
#       false = as.POSIXct(x = paste(bedtime, bt.sleep), 
#                          format = "%I:%M %p") + 
#           as.difftime(tim = 1, units = "days")), # next day for times at/after midnight
#          Wake = as.POSIXct(x = paste(waketime, wt.wake),
#                         format = "%I:%M %p") +
#           as.difftime(tim = 1, units = "days"), # time in current date plus one day, equivalently time in next day
#       `Time Difference` = difftime(time1 = Wake,
#                                    time2 = Sleep,
#                                    units = "hours")) %>%  # difference of the two times, guaranteed to be in [0, 24)
#   print(width = Inf) # TIB version 4
# 
# 
# dat.psqi$TIB <- dat.psqi$`Time Difference`
# 
# dat.psqi <- dat.psqi %>% relocate(`Time Difference`, .after = last_col())
dat.psqi <- dat.psqi %>% relocate(`TIB`, .after = last_col())
# dat.psqi <- dat.psqi %>% relocate(`Time in bed`, .after = last_col())
# dat.psqi <- dat.psqi %>% relocate(`TIB2`, .after = last_col())





# write.csv(dat.psqi, "PSQI_rawData_April2020_v1.csv", row.names = FALSE)

# currentDate <- format(Sys.time(), "%F_%H-%M")
# write_csv(dat.psqi,
#         paste("PSQI_rawData_",currentDate,".csv",sep=""),
#         na = "NA",
#         col_names = TRUE, quote_escape = "double")


# colnames(dat.psqi)

  • The PSQI data has been updated, processed, clean and scored

  • NOTE: I need to add a section that discusses the PSQI cutoff. Anything over “5” is considered significant.

4.3 INITIAL DATA SETS

Create a small dataset (dat.sm) that includes only the variables of interest.

4.4 VARIANCE CHECK

Which variable(s) in the dat.sm data has potential issues with variances?

# dat.sm.smry <- by(dat.sm, dat.sm$group.factor, summary, na.rm = TRUE)
# 
# summarise(dat.sm)

#       LOOKING AT VARIANCE AMONG VARIABLES

nearZeroVar(dat.sm,
            uniqueCut = 2,
            names = TRUE)
[1] "slpDist"

5 SAMPLE DESCRIPTIVES

5.1 Group by sex for overall sample.

#       GROUP BY SEX INFORMATION
group_by_sex <- dat1 %>%
    group_by(group.factor, sex.factor) %>%
    tally()

names(group_by_sex)[names(group_by_sex) == 'group.factor'] <- 'Group'
names(group_by_sex)[names(group_by_sex) == 'sex.factor'] <- 'Sex'

gpsex1 <- spread(group_by_sex, Sex, n) %>% 
    adorn_totals(c("row", "col"))

gpsex1 %>% 
    kable(align = "c", padding = 4, caption = "Entire Sample") %>% 
    kable_styling(c("striped", "bordered"), full_width = F) %>% 
    add_header_above(c(" ", "Sex" = 2, "")) %>%
    row_spec(0, bold = T, color = "white", background = "#461B7E" ) %>% 
    row_spec(4, bold = T, background = "lightgrey") %>% 
    column_spec(4, bold = T, background = "lightgray")

#       END     #

5.2 Group by age for overall sample.

5.3 t-test Age 2 groups

  group.factor  n
1           HC 35
2         CLBP 63
3           FM 25
  group.factor  n
1           HC 35
2         CLBP 63


    Welch Two Sample t-test

data:  age.dat$age_visit_1 by age.dat$group
t = -4.585, df = 66.101, p-value = 2.076e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.048506  -4.738601
sample estimates:
mean in group 1 mean in group 2 
       24.79230        33.18585

5.4 MEAN AGE STUFF

means.age <- dat1 %>% dplyr::select("group.factor", "age_visit_1") 

means.age  <- filter(means.age, group.factor != TRUE)

meansAge.stats <- describeBy(means.age,
                            group = "group.factor",
                            interp = FALSE,
                            skew = FALSE,
                            ranges = FALSE,
                            fast = NULL)

meansAge.stats 

 Descriptive statistics by group 
group: HC
              vars  n  mean   sd   se
group.factor*    1 35  1.00 0.00 0.00
age_visit_1      2 35 24.79 3.86 0.65
------------------------------------------------------------ 
group: CLBP
              vars  n  mean    sd   se
group.factor*    1 63  2.00  0.00 0.00
age_visit_1      2 53 33.19 12.45 1.71
------------------------------------------------------------ 
group: FM
              vars  n  mean    sd   se
group.factor*    1 25  3.00  0.00 0.00
age_visit_1      2 22 45.18 13.95 2.97
------------------------------------------------------------ 
group: TRUE
NULL
# means.stats <- filter(means.stats,  group1 !="TRUE")

# means.stats99
# 
# means.stats99 %>% kable(align = "c", caption = "test") %>% 
# kable_styling(c("striped", "bordered"), full_width = F) %>% 
#   # add_header_above(c(" ", "Sex" = 2, "")) %>%
#   row_spec(0, bold = T, color = "white", background = "#461B7E" ) %>% 
#   row_spec(4, bold = T, background = "lightgrey")# %>% 
#   # column_spec(4, bold = T, background = "lightgray")

Means and SDs for each group.

means.data <- dat1 %>% dplyr::select("group.factor",
                       "slpQual", "slpLat", 
                       "slpDur", "slpEff", 
                       "slpDayFcn", "psqi_Global",
                       "ess_total", "bdi_total", 
                       "mcgill_total") 

means.data   <- filter(means.data, group.factor != TRUE)

means.stats <- describeBy(means.data,
                            group = "group.factor",
                            interp = FALSE,
                            skew = FALSE,
                            ranges = FALSE,
                            fast = NULL)

means.stats 

 Descriptive statistics by group 
group: HC
              vars  n mean   sd   se
group.factor*    1 35 1.00 0.00 0.00
slpQual          2 35 0.77 0.65 0.11
slpLat           3 35 0.86 0.77 0.13
slpDur           4 35 0.69 0.72 0.12
slpEff           5 35 1.20 1.37 0.23
slpDayFcn        6 35 0.26 0.44 0.07
psqi_Global      7 35 5.03 2.98 0.50
ess_total        8 35 6.09 2.73 0.46
bdi_total        9 35 3.11 4.08 0.69
mcgill_total    10 35 2.14 7.35 1.24
------------------------------------------------------------ 
group: CLBP
              vars  n  mean    sd   se
group.factor*    1 63  2.00  0.00 0.00
slpQual          2 63  1.21  0.79 0.10
slpLat           3 63  1.25  0.80 0.10
slpDur           4 63  0.86  0.72 0.09
slpEff           5 63  1.35  1.31 0.16
slpDayFcn        6 63  0.22  0.42 0.05
psqi_Global      7 63  6.46  2.85 0.36
ess_total        8 63  6.84  3.50 0.44
bdi_total        9 63  9.63  7.72 0.97
mcgill_total    10 63 24.83 15.03 1.89
------------------------------------------------------------ 
group: FM
              vars  n  mean    sd   se
group.factor*    1 25  3.00  0.00 0.00
slpQual          2 25  1.76  0.72 0.14
slpLat           3 25  2.28  0.98 0.20
slpDur           4 25  1.52  0.96 0.19
slpEff           5 25  1.72  1.31 0.26
slpDayFcn        6 25  0.08  0.28 0.06
psqi_Global      7 25  9.88  3.64 0.73
ess_total        8 25  7.28  3.94 0.79
bdi_total        9 25 15.40  9.11 1.82
mcgill_total    10 25 35.16 13.60 2.72
------------------------------------------------------------ 
group: TRUE
NULL
# means.stats <- filter(means.stats,  group1 !="TRUE")

# means.stats99
# 
# means.stats99 %>% kable(align = "c", caption = "test") %>% 
# kable_styling(c("striped", "bordered"), full_width = F) %>% 
#   # add_header_above(c(" ", "Sex" = 2, "")) %>%
#   row_spec(0, bold = T, color = "white", background = "#461B7E" ) %>% 
#   row_spec(4, bold = T, background = "lightgrey")# %>% 
#   # column_spec(4, bold = T, background = "lightgray")

#PLOTS

5.5 CORR PLOTS

The end of corr plots

5.6 BAR PLOTS

5.7 BOX PLOTS

6 ANOVA (Main effect of group)

# construct anova model for each column of dat1[, my_outcome_variables]
model_summaries <- dat1[,c("age_visit_1","slpQual", "slpLat", "slpDur",
                           "slpEff", "slpDist","slpMeds",
                           "slpDayFcn", "psqi_Global","ess_total",
                           "isi_total", "bdi_total", "mcgill_total",
                           "TIB")] %>%
  purrr::map(~ aov(.x ~ group.factor, data = dat1))  %>%
  # append the TukeyHSD CI estimates
  purrr::map(function(x) {
    list(
      model = x,
      tukey = TukeyHSD(x)
    )
  })

####        THIS RETURNS THE RESULTS OF THE OVERALL ANOVA   ####
GetModels <- function(l) l$model
Models <- purrr::map(model_summaries, GetModels) %>% 
  purrr::map_dfr(broom::tidy, .id = "Name")
# head(Models)

#       
Models2 <- filter(Models,  term !="Residuals")
Models2 <- Models2 %>% 
    dplyr::select("Name" , "p.value")


####            ANOVA TABLE MADE WITH KABLE         ####
Models2 %>%
    mutate(p.value = 
            cell_spec(
                (formatC(x = p.value, 
                         digits=3, width = 3, 
                         format='g', flag = "-", drop0trailing = TRUE, 
                         preserve.width = "common")), "html",
                color = ifelse(p.value <= 0.054, 
                               yes="red", no="blue"), escape = F)
    ) %>% 
    mutate_if(is.numeric, function(x) {
        cell_spec(
            (formatC(x,digits=3, width=3, 
                     format="f", flag="-", drop0trailing = FALSE)),
            "html")}
    ) %>%
    kable("html", escape = FALSE, align = "l", padding = 4,
          caption = "ANOVA full model") %>%
    kable_styling(c("striped", "bordered", "hover"), full_width = F) %>%
    footnote(general = "Red = Significant main effect for group (p ≤ 0.05)")  %>% 
    row_spec(0, bold = T, color = "white", background = "#461B7E" )
ANOVA full model
Name p.value
age_visit_1 2.89e-09
slpQual 6.8e-06
slpLat 6.61e-09
slpDur 0.000173
slpEff 0.317
slpDist 0.287
slpMeds 2.04e-05
slpDayFcn 0.215
psqi_Global 7.29e-08
ess_total 0.374
isi_total 0.000299
bdi_total 9.01e-09
mcgill_total 1.61e-17
TIB 0.995
Note:
Red = Significant main effect for group (p ≤ 0.05) 
##      END     ##
#       USING PSQI RAW COMPONENT SCORES
# construct anova model for each column of dat1[, my_outcome_variables]
model_summary_raw <- dat1[,c("age_visit_1", "slpQualraw", "slpLatraw", "slpDurraw",
                           "slpEffraw", "slpDistraw","slpMedsraw",
                           "slpDayFcnraw", "psqi_Globalraw","ess_total",
                           "isi_total", "bdi_total", "mcgill_total",
                           "TIB")] %>%
  purrr::map(~ aov(.x ~ group.factor, data = dat1))  %>%
  # append the TukeyHSD CI estimates
  purrr::map(function(x) {
    list(
      model = x,
      tukey = TukeyHSD(x)
    )
  })

####        THIS RETURNS THE RESULTS OF THE OVERALL ANOVA   ####
# GetModels <- function(l) l$model
Models3 <- purrr::map(model_summary_raw, GetModels) %>% 
  purrr::map_dfr(broom::tidy, .id = "Name")
# head(Models)

#       
Models3 <- filter(Models3,  term !="Residuals")
Models3 <- Models3 %>% 
    dplyr::select("Name" , "p.value")


####            ANOVA TABLE MADE WITH KABLE         ####
Models3 %>%
    mutate(p.value = 
            cell_spec(
                (formatC(x = p.value, 
                         digits=3, width = 3, 
                         format='g', flag = "-", drop0trailing = TRUE, 
                         preserve.width = "common")), "html",
                color = ifelse(p.value <= 0.054, 
                               yes="red", no="blue"), escape = F)
    ) %>% 
    mutate_if(is.numeric, function(x) {
        cell_spec(
            (formatC(x,digits=3, width=3, 
                     format="f", flag="-", drop0trailing = FALSE)),
            "html")}
    ) %>%
    kable("html", escape = FALSE, align = "l", padding = 4,
          caption = "ANOVA full model with raw scores") %>%
    kable_styling(c("striped", "bordered", "hover"), full_width = F) %>%
    footnote(general = "Red = Significant main effect for group (p ≤ 0.05)")  %>% 
    row_spec(0, bold = T, color = "white", background = "#461B7E" )

##      END     ##

6.1 POSTHOC

GetTukeys <- function(m) m$tukey
Tukeys  <- purrr::map(model_summaries, GetTukeys)  %>% 
    purrr::map_dfr(broom::tidy, .id = "Name")
# head(Tukeys)
# Tukeys

names(Tukeys)[names(Tukeys) == 'Name'] <- 'Var'
names(Tukeys)[names(Tukeys) == 'estimate'] <- 'diff'



####            PostHoc TABLE USING KABLE           ####
select(Tukeys, -term, -Var) %>%
    mutate(adj.p.value = 
            cell_spec((formatC(x = adj.p.value, 
                               digits=3, width = 3, 
                               format='g', flag = "-", drop0trailing = TRUE, 
                               preserve.width = "common")), "html",
                      color = ifelse(adj.p.value <= 0.054, 
                                   yes="red", no="blue"), escape = F)) %>% 
    mutate_if(is.numeric, function(x) {
        cell_spec((formatC(x,digits=3, width=3, 
                           format="f", flag="-", drop0trailing = FALSE)),
                  "html")}) %>%
    kable("html", escape = FALSE, align = "l", 
          caption = "PostHoc group comparisons: **See** *Note* at table end",
          padding = 2) %>%
    kable_styling(c("striped", "bordered", "hover"), full_width = F) %>%
    footnote(general = "Red = Significant group tests (p ≤ 0.05)")  %>% 
    row_spec(0, bold = T, color = "white", background = "#461B7E" ) %>% 
    pack_rows(index = c(" " = 0, 
                        "age_visit_1" = 3,
                        "Sleep Quality" = 3,
                        "Sleep Latency" = 3,
                        "Sleep Duration" = 3,
                        "Sleep Efficacy" = 3,
                        "Sleep Disturbance" = 3,
                        "Sleep Medications" = 3,
                        "Sleep Daytime Functioning" = 3,
                        "PSQI Global" = 3,
                        "ESS Total" = 3, 
                        "ISI Total" = 3,
                        "BDI Total" = 3,
                        "McGill Total" = 3,
                        "Time in Bed" = 3)) #%>%
PostHoc group comparisons: See Note at table end
contrast null.value diff conf.low conf.high adj.p.value
age_visit_1
CLBP-HC 0.000 8.394 2.763 14.024 0.00169
FM-HC 0.000 20.386 13.352 27.419 1.2e-09
FM-CLBP 0.000 11.992 5.436 18.548 9.29e-05
Sleep Quality
CLBP-HC 0.000 0.435 0.067 0.803 0.0162
FM-HC 0.000 0.989 0.531 1.446 3.4e-06
FM-CLBP 0.000 0.554 0.141 0.967 0.00527
Sleep Latency
CLBP-HC 0.000 0.397 -0.020 0.814 0.0656
FM-HC 0.000 1.423 0.905 1.940 5.14e-09
FM-CLBP 0.000 1.026 0.559 1.493 2.36e-06
Sleep Duration
CLBP-HC 0.000 0.171 -0.215 0.558 0.545
FM-HC 0.000 0.834 0.355 1.314 0.0002
FM-CLBP 0.000 0.663 0.230 1.096 0.0012
Sleep Efficacy
CLBP-HC 0.000 0.149 -0.514 0.813 0.855
FM-HC 0.000 0.520 -0.304 1.344 0.296
FM-CLBP 0.000 0.371 -0.373 1.115 0.466
Sleep Disturbance
CLBP-HC 0.000 0.029 -0.016 0.074 0.292
FM-HC 0.000 0.029 -0.027 0.084 0.448
FM-CLBP 0.000 0.000 -0.050 0.050 1
Sleep Medications
CLBP-HC 0.000 0.286 -0.216 0.788 0.37
FM-HC 0.000 1.234 0.611 1.858 2.1e-05
FM-CLBP 0.000 0.949 0.386 1.511 0.000323
Sleep Daytime Functioning
CLBP-HC 0.000 -0.035 -0.236 0.166 0.911
FM-HC 0.000 -0.177 -0.427 0.073 0.216
FM-CLBP 0.000 -0.142 -0.368 0.083 0.297
PSQI Global
CLBP-HC 0.000 1.432 -0.099 2.962 0.0719
FM-HC 0.000 4.851 2.950 6.753 4.94e-08
FM-CLBP 0.000 3.420 1.703 5.136 1.85e-05
ESS Total
CLBP-HC 0.000 0.756 -0.944 2.455 0.544
FM-HC 0.000 1.194 -0.916 3.305 0.374
FM-CLBP 0.000 0.439 -1.466 2.344 0.848
ISI Total
CLBP-HC 0.000 4.333 -1.885 10.552 0.223
FM-HC 0.000 9.619 3.171 16.067 0.00198
FM-CLBP 0.000 5.286 1.668 8.903 0.00247
BDI Total
CLBP-HC 0.000 6.521 2.909 10.132 0.000109
FM-HC 0.000 12.286 7.800 16.771 5.76e-09
FM-CLBP 0.000 5.765 1.716 9.814 0.0028
McGill Total
CLBP-HC 0.000 22.683 16.179 29.186 6.59e-13
FM-HC 0.000 33.017 24.940 41.095 5.14e-14
FM-CLBP 0.000 10.335 3.043 17.626 0.00295
Time in Bed
CLBP-HC 0.000 0.228 -5.369 5.825 0.995
FM-HC 0.000 0.243 -6.709 7.194 0.996
FM-CLBP 0.000 0.015 -6.260 6.290 1
Note:
Red = Significant group tests (p ≤ 0.05) 
# scroll_box(height = "500px")
# GetTukeys <- function(m) m$tukey
Tukeys3  <- purrr::map(model_summary_raw, GetTukeys)  %>% 
    purrr::map_dfr(broom::tidy, .id = "Name")
# head(Tukeys)
# Tukeys

names(Tukeys3)[names(Tukeys3) == 'Name'] <- 'Var'
names(Tukeys3)[names(Tukeys3) == 'estimate'] <- 'diff'



####            PostHoc TABLE USING KABLE           ####
select(Tukeys3, -term, -Var) %>%
    mutate(adj.p.value = 
            cell_spec((formatC(x = adj.p.value, 
                               digits=3, width = 3, 
                               format='g', flag = "-", drop0trailing = TRUE, 
                               preserve.width = "common")), "html",
                      color = ifelse(adj.p.value <= 0.05, 
                                   yes="red", no="blue"), escape = F)) %>% 
    mutate_if(is.numeric, function(x) {
        cell_spec((formatC(x,digits=3, width=3, 
                           format="f", flag="-", drop0trailing = FALSE)),
                  "html")}) %>%
    kable("html", escape = FALSE, align = "l", 
          caption = "PostHoc (PSQI raw) group comparisons: **See** *Note* at table end",
          padding = 2) %>%
    kable_styling(c("striped", "bordered", "hover"), full_width = F) %>%
    footnote(general = "Red = Significant group tests (p ≤ 0.05)")  %>% 
    row_spec(0, bold = T, color = "white", background = "#461B7E" ) %>% 
    pack_rows(index = c(" " = 0, 
                        "age_visit_1" = 3,
                        "Sleep Quality raw" = 3, 
                        "Sleep Latency raw" = 3,
                        "Sleep Duration raw" = 3,
                        "Sleep Efficacy raw" = 3,
                        "Sleep Disturbance raw" = 3,
                        "Sleep Medications raw" = 3,
                        "Sleep Daytime Functioning raw" = 3,
                        "PSQI Global raw" = 3,
                        "ESS Total" = 3, 
                        "ISI Total" = 3,
                        "BDI Total" = 3,
                        "McGill Total" = 3,
                        "Time in Bed" = 3)) #%>% 
# scroll_box(height = "500px")

STOP HERE

7 ISI ANALYSES

This set of analyses uses a subset of subjects who have completed the ISI

Also create a dataset of subjects with ISI data (dat.isi).

7.1 ISI DATA VARIANCE

The variable(s) in the ISI data with variance issues are:

8 END - BEYOND HERE IS EXPLORATORY